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SUMMARY 


Our  aim  for  1992  was  lo  investigaic  whether  modern  quantum  nvechanital  methods  can 
be  useful  in  the  study  of  interfaces  between  dissimilar  matenals,  and  to  complete  at  least  a  mod 
est  piece  of  work  to  demonstrate  feasibility.  This  has  been  achieved. 

We  developed  an  approximate  quantum  mechanical  approach  to  the  problem  oi  dissimi 
lar  transition-metal  interfaces.  Because  of  its  approximate  nature,  we  have  devoted  a  large  effon 
to  testing  the  scheme  against  experiment  and  first-principles  theory.  Our  expertise  in  kval  den¬ 
sity  functional  theory  has  allowed  us  to  draw  upon  our  first- pnnciplcs  capability  both  in  deter 
mining  parameters  for  our  model  Hamiltonian,  and  in  testing  it. 

We  have  studied  in  some  detail  a  lattice-mismatched  interface  between  pure  Mo  and  Re 
metals  (Section  2).  Our  conclusions  concerning  this  system  are  that  the  shear  strength  of  the 
interface  is  diminished  from  the  theoretical  strength  by  a  factor  of  about  twenty;  the  reduction  is 
due  to  the  presence  of  a  misfit  dislocation  that  glides  in  the  interface  during  deformation  This 
factor  is  still  some  orders  of  magnitude  smaller  than  the  reduction  of  theoretical  strength  of  per 
feet  crystals  by  lattice  dislocations.  Therefore,  while  the  strength  of  a  mismatched  interface  is 
expected  to  be  less  than  that  of  a  lattice  matched  one  (all  other  things  being  equal),  the  misfit 
dislocation  is  not  expected  to  significantly  weaken  the  bicrystal  (Section  3). 

In  order  that  we  can  begin  to  study  systems  of  arbitrary  chemical  compositon.  such  as 
aluminide  alloys,  nonmetal  impurities  in  metal  interfaces,  and  the  metal-ceramic  interface,  we 
will  need  a  general,  first-principles,  molecular-dynamics  capability.  Part  of  our  proposal  con¬ 
cerned  our  continuing  development  of  such  a  scheme,  and  we  report  progres.s  on  this  (Section  6) 

In  conclusion,  in  the  short  time  available  we  believe  useful  results  can  emerge  from 
quantum  mechanical  calculations  that  have  some  bearing  up>on  engineering  practice.  Our  work  is 
sufficiently  complete  that  we  have  some  concrete  results  to  repiort;  at  the  same  time  there  is  much 
to  continue  in  the  same  vein  in  the  near  future. 
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AIMS  AND  ACHIEVEMENTS 


A  large  effort  is  under  way  to  understand  welded  and  diffusion-bonded  joints  between 
dissimilar  materials.  There  is  a  number  of  outstanding  questions  concerning  the  wcldabihiy  of 
traditional  engineering  materials.  Also,  newer  materials  are  emerging  and  we  want  to  learn  how 
to  bond  them  together.  For  example,  modem  ceramic  materials  allow  the  uses  of  ver)'  high 
engine  temperatures  and  we  need  to  know  how  to  join  these  to  underlying  metal  components 
Hence,  it  is  increasingly  important  to  understand  the  metal-ceramic  interface  from  a  very  funda 
mental  point  of  vicw.l  Certain  refractory  metals  arc  again  being  developed  for  applications  in 
nuclear  engines  and  as  energy  convertors  in  the  space  program,*  after  having  been  set  aside 
because  they  were  thought  less  useful  as  engineering  materials.  An  advantage  of  many  of  these 
metals  is  that  they  are  weldable.^ 

The  work  undenaken  for  the  present  ONR  contract  is  to  study  fundamental  aspects  of 
metal  joining.  The  question  is,  whether  recently  developed  quantum  mechanical  approaches  to 
materials  science  can  be  useful  in  the  context  of  the  design  of  weldable  materials.  This  is  a  high- 
risk  program  since  so  far  no-one  has  attempted  to  use  these  new  methods  in  the  context  of  such 
complicated  systems.  However,  the  payoffs  can  be  very  great  if  this  research  will  lead  to  our 
goal  of  first-principles  criteria  of  weldability. 

Certainly,  the  understanding  of  the  strengths  of  dissimilar  metal  interfaces  requires  an 
investigation  of  the  quantum  mechanical  nature  of  the  bonding  between  atoms.  A  purely  classi¬ 
cal  approach  cannot  arrive  at  a  proper  description,  since  the  chemical  bond  is  a  quantum  mechan¬ 
ical  phenomenon.  Therefore  we  require  that  our  atomistic  simulations  of  interface  structure  be 
quantum  mechanically  based,  even  if  they  cannot  be  completely  “first-principles.” 

We  intended  that  at  the  end  of  the  first  year  (1992)  we  would  have  completed  a  feasibility 
study  as  well  as  arming  ourselves  with  suitable  results  to  make  the  first- year  report  describe  a 
rather  complete  piece  of  research.  This  we  have  accomplished.  We  chose  to  pursue  the  follow  ¬ 
ing  line  of  research: 

1 .  Decide  on  a  dissimilar  metal  interface  to  study. 

2.  Develop  an  approximate  quantum  mechanical  description  (Hamiltonian)  for  Ixjth 
metals  and  their  alloys,  which  includes  a  prescription  for  interatomic  forces. 

3.  Use  molecular  statics  simulations  to  find  the  equilibrium,  zero-temperature  strength, 
and  structure  of  the  interface. 

4.  Evolve  a  plan  of  future  work  on  this  system. 

These  four  items  are  discussed  in  the  following  four  sections. 

Our  approach  has  been  to  use  empirical  quantum  theory  because  current  first-principles 
methods  are  too  computer-intensive  to  be  app.'icable  to  problems  involving  more  than  a  few  ten.s 
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of  inequivalent  atoms.  Funhermore.  the  calculations  of  interatomic  forces  are  still  problematical. 
Both  these  difficulties  can  be  removed  using  a  new  band-structure  method  that  is  currently  under 
development.  We  briefly  repon  progress  on  this  in  Section  6. 
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2. 


THE  Mo/Re  INTERFACE 


As  mentioned,  there  is  renewed  interest  in  refractory  metals  for  use  in  high- temperature 
applications  under  conditions  of  hydrogen  environments  or  neutron  fluxes.  For  example.  Mo, 

Re,  and  Mo-Re  alloys  are  strong  and  ductile  as  well  as  weldable.  We  have  chosen  to  investigate 
the  strength  of  a  Mo/Re  interface. 

This  work  takes  us  one  step  further  in  complexity  from  ideally  lattice-matched  interfaces, 
since  the  lattice  constants,  and  indeed  crystal  structure,  are  different  in  these  two  metals.  How¬ 
ever,  the  c-axis  in  hep  Re  is  4.46A  while  the  length  of  the  face  diagonal  <2(  1  lOj  in  bcc  Mo  is 
4.45A.  Therefore  an  interface  in  which  these  two  directions  arc  parallel  and  in  the  interface 
plane  will  be  lattice-matched  at  least  in  that  direction.  The  Mo(l  10)/Rc(I0l0)  interface  contains 
these  two  directions,  and  if  they  are  oriented  parallel  to  each  other,  the  orthogonal  direction  in  the 
interface  is  parallel  to  the  a-zxes  whose  lengths  are  3.15A  and  2.76A  in  Mo  and  Re.  respectively. 
The  misfit  along  this  direction  is  almost  exactly  in  the  ratic  7/8  so  that  the  interface  repeats  itself 
every  7  Mo  lattice  constants  and  every  8  Re  lattice  constants,  as  illustrated  in  Figure  1 

In  the  general  case  of  a  lattice-mismatched  interface,  the  mismatch  occurs  in  two  onhog- 
onal  directions  and  the  expected  equilibrium  structure  comprises  a  grid,  or  two-dimensional  net¬ 
work,  of  misfit  dislocations.  This  is  a  complicated  case  to  treat;  we  simplify  it  by  choosing  an 
interface  that  is  lattice-mismatched  in  only  one  direction.  The  expected  structure  is  then  a  one¬ 
dimensional  array  of  misfit  dislocations  similar  to  an  array  of  edge  dislocations  in  a  low-angle- 
tilt  grain  boundary.  We  expect  that  this  structure  will  provide  much  of  the  essential  features  of  a 
lattice-mismatched  interface  except  for  the  properties  of  dislocation  intersections. 

Our  strategy  is  as  follows.  Taking  the  atomic  arrangement  shown  in  Figure  Kb) — the 
unrelaxed,  or  epitaxial  interface — if  we  allow  atoms  to  move  according  to  the  forces  acting  on 
them,  how  will  the  structure  relax?  Presumably  the  lattice  mismatch  will  contrive  to  concentrate 
itself  at  the  core  of  a  misfit  dislocation,  which  will  be  an  edge  dislocation  lying  in  the  interface 
plane  with  its  line  sense  ^  parallel  to  the  z-direction,  normal  to  the  plane  of  the  paper  in  Fig¬ 
ure  1(b).  Having  done  this  relaxation,  or  atomistic  simulation  we  will  initially  ask  two  questions: 

1 .  What  is  the  width  of  the  core? 

2.  What  is  the  size  of  the  Peierls  barrier  for  dislocation  motion?  In  effect,  what  is  the 
shear  strength  of  this  ideal  interface?  Later,  we  can  then  ask  such  questions  such  as 
What  is  the  effect  of  interdiffusion  on  the  strength?  and  What  is  the  effect  of  the  point 
defect  dislocation  interaction? 
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(a)  Juxtapostion  of  the  Re  hep  and  Mo  bcc 
lattices.  Planes  parallel  to  the  interface 
are  shown  shaded.  The  /-direction  is  par¬ 
allel  to  [1  i  OjMo  ana  {0001  Jrb  m  this  direc¬ 
tion.  the  structure  repeats  itself  every 
repeat  of  the  hep  c-axis,  which  is  equal  in 
length  to  V2aMo-  The  x-direction  also  lies 
in  the  interlace  and  is  parallel  to  the 
a-axes  of  Mo  and  Re  which  differ  in  length 
in  the  ratio  8/7;  the  interface  is  therefore 
mismatched  in  this  direction  txjt  repeats 
every  7  Mo  and  every  8  Re  lattice  spac- 
ings  a.  The  /-direction  is  normal  to  the 
interface  and  parallel  to  [1 1  0)mo  and 
(tOiOlRo-  The  distance  n  between  the  two 
crystals  is  shown  as  a  broken  line  The 
total  energy  will  be  minimized  with  respect 
to  h. 


ooooooooooooooo 
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(b)  A  projection  of  the  Mo/Re  interface  down 
the  z-axis  of  Figure  1(a).  Only  two  planes 
of  atoms  need  be  shown  in  this  direction  as 
the  structure  then  repeats  itself.  These  two 
planes  are  shown  as  large  and  small 
atoms.  The  x-axis  points  to  the  right  and 
the  figure  shows  7  Mo  and  8  Re  lattice 
spacings  in  this  direction  after  which  the 
structure  repeats  itself.  The  /-axis  points 
upward.  Note  disregistry  in  the  interplanar 
spacing  across  the  interlace.  The  same 
Mo  atom  in  both  figures  is  slightly  shaded 
to  emphasize  the  relation  between  the  two 
perspectives. 


Figure  1 .  Mo/Re  interface. 
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3. 


TOTAL  ENERGY  AND  INTERATOMIC  FORCES 


3.1  THE  TIGHT-BINDING  FORMALISM  IN  RELATION  TO  DENSITY  FUNCTIONAL  THEORY 

As  mentioned  in  Section  2,  we  will  require  a  means  to  calculate  the  forces  acting  upon 
Mo  and  Re  atoms  in  a  solid-state  assembly  of  these  atoms.  A  traditional  approach  would  be  to 
develop  a  pair-potential  or  Finnis-Sinclaii^  model  to  describe  the  interatomic  interactions.  This 
is  not  always  possible,  particularly  as  we  shall  see  below  in  the  Mo-Re  system,  and  is  never  reli¬ 
able.  We  take  the  view  that  interatomic  interactions  are  quantum  mechanical  in  nature  and  can 
only  be  correctly  calculated  by  solving  a  Schrodinger  equation  for  the  electrons  responsible  for 
bonding  in  the  solid  state.  Our  usual  approach  is,  in  fact,  a  first-principles  approach;  we  use 
density  functional  theory  to  study  the  properties  of  transition  metals  and  alloys.  5.6  in  the  present 
situation,  however,  a  completely  first-principles  approach  is  out  of  the  question.  First  the 
ab  initio  density  functional  theory  is  too  computer  intensive  for  such  complicated  problems,  and 
second,  there  are  difficulties  involved  in  calculating  interatomic  forces  in  transition  metals. 

We  have  therefore  chosen  to  do  our  calculations  within  the  tight-binding  approximation 
to  the  electronic  structure  problem.  Our  approach  is  similar  to  that  of  Sutton  et  alP  The  way  the 
tight-binding  method  works  is  this;  We  know  which  electrons  are  most  responsible  for  bonding 
in  the  central  transition  metals  (namely  the  d-electrons),  and  we  know  the  basic  form  that  the 
Hamiltonian  must  take.  So,  instead  of  computing  the  Hamiltonian  self-consistently,  which  is 
very  time  consuming,  we  parameterize  its  matrix  elements,  after  which  we  solve  the  Schrodinger 
equation  just  as  in  the  first-principles  approach  to  obtain  eigenvalues  e/  and  eigenvectors  c,. 

In  density  functional  theory,  the  (negative)  cohesive  energy  is 

EbIP}=  fdrpj}(r)  + E^c  + Eu-Ea.om,  .  (D 

l.OCC  ^ 

The  sum  is  over  occupied  eigenvalues;  p  is  the  electron  density;  0^5  is  the  classical  electrostatic 
potential  of  the  density:* 


V2(t)ej  =  -8ll(p-p+)  ,  (2) 

where  p*  is  the  density  of  positive  charges  (nuclei);  £xc  ts  the  exchange  and  correlation  energy; 
Pxc  is  the  exchange  and  correlation  potential: 

SE^IPI 

Sp 


*  We  will  use  attomic  Rydberg  units  (a.u.)  throughout;  ttVlm  =  1  and  e^  =  2.  1  Ry  =  13.6  eV  and  the  unit  of 
length.  1  bohr  =  0.529A. 
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Eli  is  the  electrostatic  pair  repulsion  between  the  ion  cores,  and  Eatoms  is  the  total  energy  of  the 
free  atoms  in  vacuo. 

In  density  functional  theory,  e  and  c  are  eigenvalues  and  eigenvectors  of  the  effeciive, 
single-particle  Hamiltonian 

H  =  V2  +  t)^^(r) 

(in  atomic  units).  The  Schrodinger  equation, 

H  Ci  =  ei  Ci 

must  be  solved  self-consistently  with  the  Poisson  equation,  Eq.  (2),  since  the  effective  potential 


^eff  -  4^es  +  f^xc 

depends  on  the  charge  density  through  Eq.  (2). 

In  tight-binding  theory  the  cohesive  energy  usually  takes  the  general  form^-^ 

Eb  =  Ie,-^Er.p-^C  ,  (3) 

t.OCC 

where  e  are  the  eigenvalues  of  the  parameterized  tight-binding  Hamiltonian,  and  %ep  +  C  is  sup¬ 
posed  to  represent  the  rest  of  the  terms  in  Eq.  (1)  with  the  constant  C  fixing  the  zero  of  energy 
with  respect  to  the  free  atoms.  Erep  has  the  form  of  a  pair  potential.  Sutton  et  aO  were  able  to 
justify  this  form  rigorously  using  the  Harris  approximation  to  the  self-consistency  procedure.  In 
their  approach,  the  repulsive  energy  represents  the  changes  in  electrostatic  and  exchange-corre¬ 
lation  energies  in  going  from  free  atoms  to  the  solid  state  and  could  be  shown  (to  a  very  good 
approximation)  to  be  a  pairwise  repulsive  energy. 

In  problems  of  dissimilar  interfaces,  it  has  been  argued  (see  M.W.  Finnish)  that  a  com¬ 
pletely  non-self-consistent  approach  is  inappropriate.  We  therefore  use  a  self-consistent,  tight- 
binding  scheme.^  We  want  the  Hamiltonian  to  be  able  to  respond  to  changes  in  the  potential,  at 
least  in  its  diagonal  matrix  elements.  This  is  done  in  the  following  way:  After  solving  the 
Schrodinger  equation  for  the  tight-binding  Hamiltonian,  there  will  be  accumulations  of  charge 
on  each  atom  k.  These  accumulations  will  give  rise  to  an  approximate  electrostatic  potential 
(again  in  atomic  units,  =  2) 

0;  =  I  r—^Mk+2UAqj 

-n\ 

on  atom  j.  This  is  a  sum  of  two  terms: the  first  is  the  interatomic,  Madelung  potential  assum¬ 
ing  the  charge  on  each  atom  behaves  as  a  point  charge;  the  second  is  the  intraatomic  potential, 
which  has  the  opposite  sign  and  represents  the  energy  penalty  of  heaping  charge  onto  an  atom. 
The  two  terms  may  both  be  large,  but  mostly  they  tend  t~  cancel  each  other.  The  parameter 
U  is  sometimes  called  the  Hubbard  U.  is  then  the  elecv  ostatic  potential  at  atomic  sitey  that  is 
seen  by  an  electron  described  by  the  tight-binding  Hamiltonian.  We  then  adjust  the  diagonal 
elements  of  the  Hamiltonian  by  adding  the  at  each  site  and  solve  the  Schrodinger  equation 
again  to  obtain  a  new  electrostatic  potential.  This  procedure  continues  until  self-consistency  is 
reached.  In  this  way,  the  off-diagonal  matrix  elements  remain  fixed  by  their  values  determined 
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by  the  parameterization,  while  the  diagonal  matrix  elements  are  determined  self-cunsistcntiy.  At 
self-consistency,  the  electrostatic  energy  (which  must  be  included  in  the  tight-binding  energy  )  is 

^  J 

Because  of  the  large  cancellation  between  terms,  this  energy  is  usually  so  small  as  to  be  negligi¬ 
ble,  as  pointed  out  by  Harrison.  10 

Let  us  now  write  down  the  cohesive  energy  in  tight-binding  for  a  transition  metal  in 
which  we  include  only  cf-electrons  in  the  Hamiltonian.  The  off-diagonal  matrix  elements  are 
transfer  integrals  between  ^/-orbitals  on  different  sites  and  are  determined  by  parameterization,  as 
discussed  below.  The  diagonal  matrix  elements  are  determined  self-consistently.  We  also  will 
choose  a  fixed  number  of  d-electrons  for  each  transition  metal;  The  num’  er  is  for  the  atom 
type  at  the  position  labeled  k.  Electrons  will  transfer  between  sites  self-consistently  unless  a 
charge  neutrality  condition  is  imposed,^  which  amounts  to  setting  (/  =  «»;  in  reality,  U  is  closer  to 
about  10  eV.  The  input  number  of  d-electrons  will  be  used  in  setting  the  constant  C  in  Eq.  (3); 
this  is  only  constant,  however,  in  non-selfconsistent  tight-binding  where  the  diagonal  elements  of 
the  tight-binding  Hamiltonian  are  constants.^  The  cohesive  energy  is 

I  1  Wa  +  +  4  •  (4) 

i.occ  k 

Except  for  the  last  term,  this  is  the  same  expression  given  by  Paxton  and  Sutton, 

It  is  not  difficult  to  differentiate  this  expression  to  obtain  interatomic  forces.  Chadi^^  and 
Kohyama  et  give  details.  In  the  present  work,  we  do  make  the  approximation  that  only  the 
first  term  and  %ep  contribute  to  the  force  on  an  atom, 

3.2  PARAMETERIZATION  OF  THE  HAMILTONIAN 

We  need  to  find  parameters  for  the  Hamiltonian  matrix  elements  and  the  pair  potentials 
that  make  up  E^gp.  We  also  need  values  for  the  Hubbard  U  parameter.  It  is  in  this  sense  that  the 
tight-binding  method  is  empirical  rather  than  first-principles  and,  for  thi.'  ’■eason,  we  will  have  to 
test  the  parameters  we  use  to  ensure  that  they  give  sensible  results  for  the  known  properties  of 
Mo  and  Re.  We  can  then  apply  the  parameterized  Hamiltonian  to  the  problem  of  the  strength 
and  structure  of  interfaces  with  confidence. 

The  tight-binding  theory  of  transition  metals  is  well  developed  and  goes  back  some 
20  years.i5  In  fact  the  method  has  been  practically  applied  with  great  success  long  before  Sutton 
et  al?  and  Foulkes  and  Haydock^^  made  the  connection  between  tight-binding  and  density  func¬ 
tional  theories.  There  are  two  uncertainties  in  the  method.  First  is  the  assumption  that  the 
Hamiltonian  can  be  parameterized  and  that  the  remaining  energy  can  be  cast  in  pair-wise  form. 
Second  is  the  need  to  find  functional  forms  and  parameters  for  the  Hamiltonian  matrix  elements 
and  the  pair  potential.  Specifically,  these  forms  and  parameters  need  to  depend  only  on  the  dis¬ 
tance  between  the  two  atoms  in  question  and  the  types  of  atoms  involved.  Implicit  in  this  is  the 
further  assumption  that  the  Hamiltonian  is  “two-center”;  that  is,  it  has  no  terms  involving  three 
atoms  simultaneously.  This  is  closely  connected  with  a  final  assumption,  namely  that  the 
Hamiltonian  appears  in  an  orthogonal  representation.  Both  nonorthogonality  and  three-center 
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terms  are  known  to  vanish  to  first  order  m  the  band  theor)  ot  transition  mctab  Gi  ven  iluj;  irsc 
first  assumpaon  is  ngorousiy  defcnsibic  t'rom  within  density  funciionai  thcors ,  wc  arc  Ictt  wiiM 
the  selection  ot  parameters  that  will  yield  an  adequate  descnpiu.m  ot  the  transuuin  inculs  in 
ques-tion.  Nevertheless,  because  this  is  an  empirical  method,  we  need  to  make  csii.susU’.e  tesr 
of  the  model  to  show  that  u  is  suitable  for  use  m  atomistic  vitnuUmon 

Two  functional  fomis  arc  popular  for  the  Haimltonian  matrix  clcmcrns  and  tJie  pan 
potential.  In  the  canonical  mxxlcl.  the  disia.ncc  dependence  is  the  inverse  fittl!  p^swci  lur  the 
matrix  dements  and  some  larger  power  (around  lOi  fur  the  pair  fxjtcnual  While  tins  is  she 
model  that  denves  naturally  from  band  thcor>',‘ '  wc  have  found  in  the  present  work  tliat  is  taiis 
to  reproduce  structural  energy  ordenng  as  desenbed  below .  Instead,  wc  have  found  the  nuxic!  le 
Spanjaard  and  Desjonquferes*^  to  be  the  most  robust  and  reliable  In  this  mosiel,  the  twi*  center 
matrix  element  of  the  Hamiltonian  connecting  two  atoms  separated  by  a  distance  r  takes  ttie  form 
h{r)  =  and  a  pair  potential  of  the  fonn  <p<r^  -  be  P'  ts  used  Spanjaard  and  Desjonquerrs 

find  p/q  =  2.95  by  compan.son  with  the  Rose  equation  of  state  and  tabulate  values  of/?  and  q  iot 
all  the  transition  metals.  Wc  do  not  argue  that  there  is  any  physical  basis  for  their  choice  of ;? 
and  q,  but  we  do  find  that  these  are  sensible  parameters.  The  ouicoitk  of  any  calculation  must 
anyway  be  insensitive  to  details  of  the  parameter  choice,  and  this  is  index'd  the  case  Giscn  p  and 
q  for  Mo  and  Re,  wc  have  only  two  free  parameters  for  each  metal,  namely  the  prefactors  f  and  h 
We  will  adjust  b  to  obtain  the  correct  lattice  constant  (atomic  volume}  and  /  will  be  chosen  to 
obtain  the  same  width  of  the  d-band  a.s  found  in  a  sclf-consistcnt  hand  caiculaiion  f  inalis ,  the 
ratio  of  the  ddo:ddn  ddS\s  given  the  usual  canonical  fonn,  6:  4, 1 , 

We  should  point  out  that  our  mode!  includes  only  d  elcctrons  tn  the  Hamiltonian  The 
effect  of  the  5-electrons  is  to  provide  a  positive  pressure,  which  is  included  in  the  repulsive 
pairwise  term.  The  number  of  d-clectrons  S'a  is  therefore  a  parameter,  not  necccssanly  the  num 
her  of  c  electrons  in  the  free  atom.  The  cohesive  energy  in  Eq.  (4  )  is  therefore  not  determined 
strictly  with  respect  to  free  atoms  but  to  some  vacuum  reference  state  in  w  hich  the  number  of 
d-electrons  is  Nd.  The  most  sensible  choice  of  S’d  should  come  from  exact  band  calculations,  in 
which  case,  we  find  4.39  d-electrons  in  Mo  and  5. 12  in  Re  The  difference  is  less  than  one 
because  the  5d  scries  have  more  s  electrons,  in  general,  than  the  4<i  (Figure  2).  What  is  staniing 
is  that  in  Mo  the  bcc  structure  is  strongly  favored,  while  the  bcc  structure  is  very  unfasoraMc  in 
Re.  Figure  2  shows  hep-fee  and  bcc-fcc  energy  differences  calculated  using  band  theory  It  is 
well  known  that  what  determines  the  stability  of  transition  metals  is  simply  the  number  of 
d  electrons.  However,  the  rapid  ri.se  in  the  bcc-fcc  energy  difference  with  a  change  in  A  j  of  less 
than  one  cannot  be  reproduced  in  a  simple  tight-binding  mcxiel  without  5<i-hybndization  Wc 
will  return  to  some  of  the.se  points  later,  but  for  now-  we  mention  that  we  have  tned  adjusting  S'a 
and  found  that  the  best  value  for  Mo  is  the  band  theory  value  of  4.39;  however,  if  we  choose 
5.12  in  Re.  then  the  bcc  structure  is  stabilized,  so  that  to  gu.irantce  the  correct  ordenng  of  cry  si.il 
structures  we  have  to  set  iV*'  =  6. 

Finally,  we  need  to  determine  the  Hubbard  U  parameter,  and  fur  this  also  we  turn  to 
density  functional  theory.  We  use  the  LMTO-.ASA  method  in  the  following  way.  First,  we 
make  a  self-consistent  calculation  of  the  atomic  sphere  charge  density  and  logarithmic  denva- 
tives  in  Mo  and  Re  at  equlibrium  volume.  Keeping  the  boundary  conditions  fixed,  we  add  a 
small  charge  dqa  into  the  d-electron  channel  and  recompute  the  total  energy  of  the  atomic  sphere 
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Figure  2,  Energy  differences  between  the  bcc  and  fee  lattices  (squares) 
and  hep  and  fee  lattices  (hexagons)  in  the  4d  (open  symbols) 
and  5d  (closed  symbols)  central  transition  metals  (groups  5- 
9)3^  Both  the  energies  and  cFelectron  numbers  are 
determined  from  band  theory.  Note  the  rapid  destabilization  of 
the  bcc  phase  over  a  very  small  change  in  N(j  This  effect 
cannot  be  reproduced  in  a  simple  <#-band  tight-binding  model 

We  then  use 


to  determine  the  Coulomb  integral  U  for  the  d  electrons.  This  is  illustrated  in  Figure  3. 


The  parameters  for  Mo  and  Re  thus  determined  are  displayed  in  Table  1 .  (Those  for  Re 
have  been  slightly  adjusted  to  make  the  conditions  cne  =  and  a}ji^aRe  -  8/7  exact.)  It 

only  remains  to  specify  the  interactions  between  dissimilar  atoms.  For  these,  we  use  the 
geometric  means  of  the  pure  atom  interactions,  which  is  the  correct  combination  according  to 
band  theory. 


Table  1 

PARAMETERS  FOB  THE  TIGHT-BINDING  MODEL  FOR  Mo  AND  Re 
(All  numbers  are  In  atomic  units) 


P 

<7 

f 

b 

Na 

U 

Mo 

1.9513 

0.6621 

0.4679 

1408.8 

4.39 

0.854 

Re 

2.0316 

0.5900 

1893.9 

6.00 

0.872 

Figure  3.  Total  energy  in  the  local  density  approximation 
for  Mo  and  Re  as  a  function  of  a  small  excess 
rf-electron  charge.  The  curvature  of  these  curves 
for  the  neutral  solid  is  the  Hubbard  L/ parameter. 


3.3  TESTS  OF  THE  HAMILTONIAN 

As  mentioned,  we  want  to  test  the  tight-binding  model  extensively  before  applying  it  to 
the  Mo/Re  interface.  The  tests  fall  into  two  categories:  tests  of  the  cohesive  and  elastic  proper¬ 
ties  of  the  equilibrium  structure,  and  tests  of  structural  stability  in  which  the  atomic  coordination 
is  radically  different  from  the  equilibrium  structure.  We  will  also  test  the  finite  strain- 
deformation  behavior.  These  are  described  in  the  following  three  subsections. 

3.3.1  Equilibrium  and  Cohesive  Properties 

We  do  not  expect  to  be  able  to  reproduce  the  experimentalcohesive  energy  in  our  model. 
This  is  contrary  to  previous  claims'^’21  that  have  been  made.  Indeed,  the  cohesive  energy  does 
come  out  in  quite  good  agreement  with  experiment  if  one  ignores  the  spin  polarization  of  the  free 
atom.  However,  this  spin  polarization  is  very  large  (about  0.2  Ry)  in  the  central  transition  metals 
and  cannot  be  neglected.  If  it  is  included,  then  we  find  a  cohesive  energy  of  about  2/3  the  exper¬ 
imental  value.  The  remaining  1/3  is  due  to  sd  hybridization,  which  is  not  in  the  model  (see  Fig¬ 
ure  4).  Obviously,  one  could  then  ask,  why  not  include  the  i-band?  And  the  answer  is  twofold. 
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log2  s  (a.u.) 


Figure  4.  Band  theory  calculations  of  the  contributions  to  the  cohesive  energy  (top 
panel)  and  pressure  (lower  panel)  as  non-spin-polarized  Mo  atoms  are 
brought  from  infinity  in  the  vacuum  together  to  form  the  bcc  crystal.  The 
equilibrium  volume  is  indicated  by  a  vertical  dotted  line;  and  there  we  see 
that  the  cohesion  is  a  balance  of  a  band  contribution  which  is  the  first 
term  in  Eq.  (1)  and  renormalization  contribution  AEot  which  represents  the 
other  terms  in  Eq.  (1).  In  the  lower  panel,  we  see  clearly  that  the  cf-electron 
contribution  to  the  pressure  is  negative  (attractive)  and  the  contribution  from 
the  free-electron-like  ^electrons  is  positive.  As  noted  in  the  text,  this  latter 
effect  cannot  be  reproduced  in  a  simple  tight-binding  band  theory  and  is 
relegated  to  the  repulsive  pair-potential.  The  atomic  volume  is  represented 
as  the  logarithm  of  the  Wigner-Seitz  radius  s  in  this  plot;  the  contributions  to 
the  cohesive  energy  from  sp-  and  d-electrons  separate  out  simply  as  the 
areas  under  the  curves  from  equilibrium  to  infinite  volume.  For  further 
discussions,  see  Pettifor^s,  Gellot  etal.  and  Anderson  etal.^^ 
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Most  important,  the  i-band  is  free-eiectron-like  and  cannot  be  described  in  tight-bjnding  As  the 
crystal  is  compressed,  the  ^-electrons  are  forced  closer  to  the  core  where  ihetr  kinetic  energy  is 
raised  by  orthogonality  constraints.  The  bottom  of  the  5-band  then  increases  w  ith  decreasing  vol¬ 
ume  and  the  5-pressure  is  positive.  The  d-pressure  is  negative;  indeed,  in  a  simple,  orthogonal 
tight-binding  model,  one  cannot  reproduce  a  positive  band  pressure,  so  that  the  effect  of 
5-eIectrons  is  best  relegated  to  the  repulsive  pair  term.  A  second  reason  to  omit  the  5-band  is  that 
one  cannot  build  a  free-electron-like  5-band  out  of  tight-binding  5-functions.  One  would  need  to 
include  p-orbitals  as  well,  which  would  make  the  basis  too  large  for  practical  purposes. 

We  do,  however,  expect  to  reproduce  other  equilibrium  properties — for  example  elastic 
constants  and  heats  of  formation  of  alloys — in  which  the  ffee-atom  energies  cancel  out.  The 
cancellation  is  a  test  of  the  interaction  parameters  between  Mo  and  Re;  we  have  imagined  a 
hypothetical  MoRe  alloy  in  the  ordered  bcc  (CsCl  or  B2)  structure  and  used  full-potential,  self 
consistent  LMTO  theory  to  compute  its  lattice  constant,  elastic  constants,  and  heat  of  formation. 
We  have  compared  these  with  the  predictions  of  the  self-consistent  tight-binding  model.  We 
have  also  used  the  model  to  calculate  the  elastic  contants  of  bcc  Mo  and  the  axial  ratio  of  hep  Re 
and  compared  these  with  experiment.  The  results  are  displayed  in  Table  2. 

Table  2 

CALCULATED  EQUILIBRIUM  PROPERTIES  COMPARED  WITH  EXPERIMENT  OR 

EXACT  THEORY  (In  parentheses) 


(The  lattice  constant  a  is  in  bohr,  c/a  is  the  hep  axial  ratio,  C'  is  (Ci2-C44)/2  and  K  is 
the  bulk  modulus— -these  are  given  in  GPa;  Hf,  the  heat  ot  formation,  is  in 
mRy/atom.) 


a 

C/a 

C 

C44 

K 

Hf 

Mo 

5.95 

112 

126 

299 

(151) 

(109) 

(262) 

Re 

5.22 

1.610 

(1.614) 

MoRe 

5.89 

91 

143 

340 

+1.4 

(5.89) 

(89) 

(217) 

(299) 

(+8.5) 

The  agreement  is  about  what  one  would  expect  from  a  tight-binding  theory,  and  quite 
acceptable.  The  most  obvious  disrepancies  are  in  MoRe,  where  tight-binding  underestimates  C44 
and  the  heat  of  formation.  It  is  noteworthy  that  H{  is  positive;  i.e.,  the  compound  does  not  exist, 
and  this  at  least  is  correctly  reproduced.  We  note  also  that  the  elastic  constants  of  pure  transition 
metals  are  also  determined  by  subtle  band-structure  effects  that  are  not  perfectly  reproduced  by 
the  tight-binding  model. 

3.3.2  Structural  Stability 

It  is  of  extreme  importance  to  test  the  stability  of  the  tight-binding  model  against  the 
stability  of  competing  structures,  both  close-packed  fee,  hep,  and  bcc  and  also  non-clo.se-packed 
structures  such  as  simple  hexagonal  (hex)  and  simple  cubic  (sc). 5-20  Energy-volume  curves  for 
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these  structures  for  the  }d  transition  series  are  available^  for  comparison,  [•igurc  5  shou  s  ihe 
curves  calculated  using  the  tight-binding  mtxiel.  Apart  from  absolute  values  of  the  ci>hcsne 
energies,  which  are  underestimated  (as  discussed  earlier)  owing  to  lack  of  .vj-hybndi/ation  in  the 
model,  the  curves  agree  with  first-principles  theory'  very  well.  In  panicular,  the  energy  differ 
ences  between  close-packed  phases  are  in  excellent  agreement  with  band  theory'  (Figure  2i:  the 
comparison  is  shown  in  Table  3. 

Tables 

ENERGY  DIFFERENCES  (In  mRy/atom) 

BETWEEN  THE  CLOSE-PACKED  PHASES  OF 
Mo  AND  Re 


(Results  from  the  fight-binding  model  are  compared 
with  first-principles  band  theory,’®  in  parentheses  ) 


bcc-fcc 

hep-fee 

Mo 

-27.1 

+3  0 

(-34.0) 

(+2.6) 

Re 

■t-28  3 

-2  2 

(+26.4) 

(-5.9) 

Figure  5.  Cohesive  energy  versus  atomic  volume  Q.  Qq  'S  the  equilibrium 
volume.  The  curves  for  Mo  and  Re  are  computed  using  the  tight- 
binding  model.  Apart  from  errors  in  the  absolute  cohesive 
energy  due  to  neglect  of  sd-hybrkJization,  they  are  in  very  good 
agreement  with  what  one  would  expect  from  a  first-principles 
calculation.^ 
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3.3.3  Finite  Deformations 


As  a  final  test  of  our  model,  we  calculate  the  energy  as  a  function  of  homogeneous  lattice 
deformation.  Here,  the  bcc  lattice  is  given  a  simple  shear  on  the  (112)  plane  in  the  1 1 1 1 J  direc¬ 
tion  until  the  lattice  shears  into  itself.  This  is  a  homogeneous  tw'inning  shear.  Firsi-pnnciple.s 
calculations  of  the  energy  versus  shear  are  available  for  comparison,^^  which  is  shown  m  Fig¬ 
ure  6. 


Figures.  Energy  density,  W,  versus  shear  x/s  For  details  see  the 
first-principles  calculations  of  Paxton  el  al.^.  The 
agreement  between  tight-binding  and  first-principles 
calculations  is  excellent  in  the  twinning  branch,  but  rather 
poorer  in  the  antitwinning  branch. 

At  last,  we  find  quite  a  large  discrepancy  between  the  model  and  first  principles,  namely 
the  antitwinning  branch  of  the  energy-shear  curve  in  Mo.  This  is  a  very  large  distortion  of  the 
lattice  away  from  equilibrium,  so  the  disagreement  is  not  too  disturbing.  It  is  well,  however,  that 
we  have  pushed  our  model  far  enough  to  find  the  point  at  which  it  eventually  breaks  down.  Even 
classical  potentials  are  rarely  tested  as  thoroughly  as  we  have  tested  the  tight-binding  model  here. 
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4. 


RELAXATION  AND  ATOMIC  STRUCTURE 


Armed  with  a  quantum  mechanical  model  of  the  Mo-Re  system,  we  go  ahead  and  com¬ 
pute  the  atomic  structure  of  our  chosen  interface.  This  is  done  by  setting  up  the  unrelaxcd  inter- 
face  in  [Figure  1(b)],  calculating  the  forces  acting  between  the  atoms,  and  relaxing  their  positions 
in  a  molecular  statics  algorithm. '2  Even  using  a  simplified  quantum  mechanical  model,  this  is  a 
very  time-consuming  business. 

We  do  our  calculations  in  supercells  so  that  we  can  exploit  Bloch’s  theorem  m  evaluating 
the  eigenvalue  sum  in  Eq.  (4).  We  do  this  by  a  tetrahedron  integration  over  the  Brillouin  zone, 
using  72  ^-points  in  the  whole  zone.  The  size  of  the  supercells  is  such  that  there  arc  eight  layers 
of  Re  atoms  and  five  layers  of  Mo  before  the  cell  repeats  itself  along  the  y-direction  (Figure  1 ) 
We  have  done  some  tests  using  larger  cells  and  convinced  ourselves  that  the  electronic  structure 
is  converged  at  this  supercell  size.  Artificial  elastic  interactions  will  remain  and  persist  to  very' 
large  cell  sizes,  and  we  have  to  live  with  them. 

Because  of  the  complexity  of  the  problem,  we  arc  forced  to  make  some  further  approxi¬ 
mations.  The  separation  h  between  the  crystals  should,  ideally,  be  relaxed  along  with  the  atom 
positions.  Instead,  we  have  determined  h  for  the  unrelaxed  cell  and  used  that  value  (4.63  a.u.)  as 
fixed.  All  relaxations  subsequently  are  at  constant  volume.  We  also  do  the  self-consistency  pro¬ 
cedure  separately  from  the  atomic  relaxations.  That  is  to  say,  we  first  relax  the  positions  and 
then  make  the  Hamiltonian  self-consistent;  we  have  found  that  further  relaxations  are  not  signifi¬ 
cant.  Finally,  we  calculate  the  total  energy  of  the  supercell  using  twice  the  number  of  /t-points. 

The  range  of  the  Hamiltonian  matrix  elements  A(r)  and  the  pair  potential  <p(r)  are  usually 
taken  as  rather  short  in  tight-binding  theory.  However,  we  have  chosen  not  to  truncate  these  until 
after  fifth  neighbors  in  the  bcc  structure  (at  -v/Sabcc)  and  sixth  neighbors  in  the  hep  structure  (at 
2ahcp)-  This  actually  improved  the  band  structure  somewhat;  more  important,  it  devolves  from 
us  the  responsibility  of  choosing  what  arc  and  are  not  neighbors  since,  at  the  cut-off,  the  matrix 
elements  are  insignificantly  small.  Because  we  are  using  a  k-space  method,  the  longer  range  of 
the  interactions  does  not  affect  the  computing  time;  in  a  real-space  approach  such  as  the  recur¬ 
sion  method  this  would  not  be  economical. 

4.1  THE  MISFIT  DISLOCATION 

The  relaxed  structure  of  the  Mo/Re  interface  from  Figure  1(b)  is  shown  in  Figure  7.  As 
expected,  to  the  left  and  right  pans  of  the  interface,  the  lack  of  registry  (disregistry)  between 
atomic  planes  across  the  interface  has  improved  at  the  expense  of  the  appearance  of  a  dislocation 
core  at  the  center.  In  other  words,  the  disregistry,  which  is  evenly  distributed  across  the  epitaxial 
interface,  has  become  localized  at  the  core  of  a  misfit  dislocation  after  atomic  relaxation.  As  far 
as  we  know,  this  is  the  first  quantum  mechanical  simulation  of  a  dislocation  core  structure  in 
transition  metals.  We  therefore  regard  this  as  a  significant  step  forward.  Already  we  can  answer 
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Figure  7.  Structure  of  the  relaxed  interface.  The  shaded 
atom  corresponds  to  that  in  Figure  1(b)  and 
appears  at  the  core  of  the  misfit  dislocation. 

our  first  question  of  Section  2,  as  to  the  width  of  the  core.  The  width  is  roughly  half  the  disloca¬ 
tion  spacing  (there  is,  of  course,  an  infinite  array  of  dislocations  generated  by  the  periodicity 
along  the  j:-axis). 

Let  us  make  a  few  comments  about  Figure  7.  The  relaxation  of  the  Mo  layer  is  confined 
to  some  bending  of  the  atomic  planes  about  the  z-axis  due  to  long-range  elastic  effects.  These 
are  artifacts  of  the  supercell  size  mentioned  above.  The  distortions  at  the  core  are  otherwise  con¬ 
fined  to  the  Re  layer  which  moves  closer  to  the  Mo  in  the  region  of  good  registry,  and  away  from 
the  Mo  at  the  dislocation  core.  Therefore,  contrary  to  the  Peicrls  model,  there  is  some  bending  of 
the  atomic  planes. 

The  position  at  the  interface  indicated  with  a  shaded  Mo  atom  in  Figure  1(a)  turns  out  to 
be  the  center  of  the  core.  This  is  in  fact  the  lowest-cncrgy  position  for  the  dislocation,  as  we 
shall  see  below,  indicating  the  position  of  the  Peierls  valley.  In  this  alignment.  Mo  and  Re  atoms 
face  each  other  across  the  interface  in  the  same  atomic  plane  along  the  z-axis.  The  interface  pre¬ 
fers  to  remain  in  registry  at  positions  where  Mo  and  Re  atoms  face  each  other  in  atomic  planes  at 
different  distances  along  the  z-axis.  In  other  words,  referring  to  Figure  1  (b),  the  interface  prefers 
to  maintain  the  middle  region  in  registry  and  concentrate  the  misfit  at  the  edge  regions.  Put  yet 
another  way,  referring  to  Figure  1(a)  the  body-centering  atom,  rather  than  sitting  opposite  the 
center  of  the  (1010)  face  (as  shown),  prefers  to  sit  along  the  edge  of  that  face.  This  is  because  in 
the  situation  shown  in  Figure  1(a),  the  shaded  Mo  atom  is  closer  to  its  opposite  Re  atom  than  the 
nearest-neighbor  distance  in  the  pure  metals.  The  edges  of  the  interface  in  Figure  1(b)  are  there¬ 
fore  the  high-energy  regions  and  the  center  is  a  low-energy  region  since  opposite  Mo  and  Re 
atoms  are  in  different  layers  along  the  z-axis  and  therefore  not  too  close  together.  Thus,  when 
the  structure  relaxes,  the  high-energy  part  becomes  the  core  of  the  dislocation  while  the  the  low- 
energy  part  remain  in  registry.  Figure  1(b)  also  represents  the  lowest-energy  configuration  for 
the  dislocation.  The  core  lies  in  a  Peierls  valley,  since  if  it  is  to  move  in  the  x-direction  the 
disregistry  must  shift  into  the  region  of  low  energy,  which  costs  energy.  Section  4.2  discusses 
what  happens  when  we  try  to  move  the  dislocation. 
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4.2  THE  PEIERLS  ENERGY 

In  a  series  of  simulations  in  molecular  statics,  we  apply  successive  shear  strains  to  the 
atomic  cell.  It  is  clear  from  a  study  of  Figure  1(b)  that  if  we  displace  the  (upper)  Re  bkx;k  with 
respect  to  the  (lower)  Mo  block  by  a  vector  of  length  b  =  umo  -  =  0-73  bohr  in  the 

;t-direction,  then  we  would  recover  an  identical  atomic  structure.  This  is  the  translational  peri¬ 
odicity  of  the  interface  in  the  x-direction.  In  the  process,  the  dislocation  will  move  through  a 
distance  b,  achieving  a  plastic  strain  of  the  interface.  The  energy  barrier  involved  in  this  process 
is  the  Peierls  energy  of  the  misfit  dislocation. 

Beginning  with  the  relaxed  smicturt  shown  in  Figure  7,  we  apply  a  strain  of  0.2h  by 
adjusting  the  boundary  conditions  of  the  supercell  (see  Paxion26  for  details  of  how  this  is  done). 
We  again  relax  the  atomic  positions,  make  the  Hamiltonian  self-consistent  and,  obtain  the  total 
energy  of  the  cell  for  this  displacement.  Using  the  new,  sheared  structure,  we  apply  another 
strain  of  0.1  h  and  repeat  the  procedure.  We  continue  to  shear  the  cell  by  increments  of  0.  id, 
relaxing  the  positions  and  the  Hamiltonian  at  each  step.  We  obtain  a  plot  of  the  energy  of  the 
cell  versus  shear,  as  shown  in  Figure  8.  We  show  displacements  up  to  only  0.4h;  for  larger  val¬ 
ues  we  have  encountered  problems  in  determining  the  global  minimum  from  a  number  of 
metastable  states.  By  symmetry,  however,  we  know  that  the  curve  must  be  symmetrical  about 
bf2  so  only  that  value  of  the  displacement  is  lacking  in  the  present  work. 

We  can  arrive  at  a  very  important  conclusion  from  these  calculations.  As  the  dislocation 
moves  from  one  Peierls  valley  to  the  next  (a  distance  b  away)  the  energy  barrier  A£  is  on  the 
order  of  3  mRv  The  displacement  is  bf2  so  that  the  stress  is,  very  approximately, 

2A£ 

where  Aceli  is  the  area  of  the  interface  in  the  unit  cell.  This  is  about  0.4  GPa  which  is  more  than 
10  times  smaller  than  the  theoretical  strength,  which  must  be  on  the  order  of  10-20  GPa.  There¬ 
fore,  while  the  presence  of  the  dislocation  clearly  reduces  the  strength  of  the  interface,  it  does  not 
have  as  strong  an  effect  as  a  lattice  dislocation,  which  is  expected  to  reduce  the  strength  of  a 
crystal  by  several  orders  of  magnitude  from  the  theoretical  strength. 

Because  the  tight-binding  model  can  reproduce  the  theoretical  strength  of  Mo  very  well 
compared  with  a  first-principles  calculation  (see  Figure  6)  we  regard  the  above  result  as  signifi¬ 
cant.  We  repeat  our  conclusion:  the  misfit  dislocation  has  the  effect  of  reducing  the  strength  of 
the  interface  from  the  Frenkel  theoretical  strength,  although  not  by  as  much  as  would  be  expected 
from  a  lattice  dislocation  in  a  single  crystal.  Therefore,  the  presence  of  misfit  dislocations  in  an 
otherwise  perfect  interface  is  not  expected  to  reduce  the  strength  significantly.  On  the  other  hand 
there  is  a  notable  reduction  in  the  theoretical  strength  compared  with  a  lattice  matched  interface. 
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Figure  8.  Energy  f  of  the  interface  relative  to  that  of  the 
undisplaced  structure  E^of  Figure  7.  The 
displacement  is  in  units  of  the  translational  period 
b  of  the  interface. 
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5. 


PLAN  OF  FUTURE  WORK 


We  are  confident  that  the  simple  tight-binding  model  is  adequate  for  the  study  of  central 
transition-metal  interfaces.  If  we  want  to  study  transition-metal  aluminide  alloys,  metal-ceramic 
interfaces,  or  nonmetal  impurities,  then  we  would  need  to  use  first-principles  density  functional 
theory.  As  we  will  see  (Section  6)  our  program  of  work  to  develop  a  first-principles  molecular 
dynamics  for  all  elements  in  the  Periodic  Table  is  sufficiently  far  ahead  for  us  to  regard  this 
extension  as  becoming  feasible  in  the  coming  year. 

Meanwhile,  there  is  plenty  we  can  learn  about  dissimilar  metal  interfaces  in  general  from 
our  tight-binding  work.  Our  plan  for  future  work  includes  the  following: 

1.  We  would  like  to  extend  our  implementation  of  the  tight-binding  method  to  include 
molecular  dynamics  as  well  as  statics.  This  can  be  easily  done. 

2.  Now  that  we  know  the  structure  and  have  an  estimate  of  the  strength  of  the  ideal 
Mo/Re  interface,  we  need  to  investigate  the  effect  of  interdiffusion.  This  is  the  real¬ 
istic  case  encountered  in  diffusion  bonding,  and  we  want  to  know  how  the  intermix¬ 
ing  of  the  atoms  affects  the  properties  of  the  interface. 

3.  We  have  a  way  of  studying  interface  dislocations  in  general,  now  that  we  have  this 
nice  model  of  a  dislocation  which,  since  it  is  an  interfacial  dislocation,  has  no  long- 
range  strain  field.  We  would  like  to  use  this  model  system  to  study  the  dislocation 
point-defect  interaction  which  has  not  yet  been  investigated  using  a  serious  quantum 
mechanical  model. 
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6. 


ADVANCES  TOWARDS  REAL-SPACE,  FIRST-PRINCIPLES 
MOLECULAR  DYNAMICS 


As  we  have  mentioned,  the  tight-binding  method  is  suitable  only  in  certain  favorable 
cases,  for  example  silicon,  germanium,  and  the  central  transition  metals.  Eventually  we  want  to 
be  free  to  study  any  system,  including  grain  boundary  impurities  and  the  metal-ceramic  interface. 
To  this  end,  some  of  our  time  under  this  contract  has  been  spent  on  continuing  development  of 
our  new  real-space,  total-energy  method,  which  is  mosdy  being  supported  under  a  separate  ONR 
contract.  In  this  final  section,  we  repon  on  our  progress.  We  begin  by  reviewing  the  state  of  the 
art  in  first-principles  molecular  dynamics. 

Extensive  research  has  now  established  that  it  is  possible  to  calculate  in  density  func¬ 
tional  theory  (DFT),  and  with  fairly  high  accuracy  (often  comparable  to  tiic  uncertainties  in  mea¬ 
surement),  a  broad  range  of  mechanical  and  structural  properties  of  most  solids.  DFT  is,  in 
principle,  exact:  solving  it  is  completely  equivalent  to  solving  the  many-body  Schrbdinger 
equation.  The  exact  functional  is  unknown,  but  the  important  advances  in  this  area  have  been 
made  using  a  local  approximation  to  the  exact  functional.  The  local  density  approximation 
(LDA)  contains  no  empirical  or  adjustable  parameters,  but  makes  it  possible  to  obtain  an  explicit 
form  for  the  Schrbdinger  equation  that  recasts  the  problem  as  a  set  of  independenf  electrons  mov¬ 
ing  in  an  effective  one-electron  potential.  Its  success  in  predicting — from  first  principles — 
mechanical  properties  such  as  structural  energy  differences,  defect  energies,  and  shear  moduli  is 
remarkable. 

Solution  of  the  LDA  is  technically  difficult  and  computationally  demanding.  While 
many  problems  have  been  successfully  broached  within  the  LDA,  in  practice  limits  to  computa¬ 
tional  capacity  sharply  circumscribe  what  it  can  feasibly  address.  Tl.i'  difficulty,  one  of  practice 
rather  than  of  principle,  is  compounded  because  the  computational  effort  grows  rapidly  with  the 
complexity  of  the  calculationj^  The  one-electron  Hamiltonian  H  is  projected  orto  a  basis  set.  and 
obtaining  the  eigenvalues  of  H  amounts  to  the  solution  of  a  linear  algebraic  eigenvalue  problem. 
The  size  of  the  basis,  and  thus  of  H,  grows  in  proportion  to  the  number  of  inequivalent  atoms  n. 
making  the  operations  count  increase  as  the  cube  of  n.  It  is  now  routine  to  study  systems  with 
fewer  than  twenty  or  so  inequivalent  atoms,  but  the  barrier  rapidly  sets  in  beyond  this  point. 
Even  this  number  is  quite  sufficient  to  study  in  detail  many  materials  properties,  such  as  those 
needed  to  compute  phase  diagrams.  On  the  other  hand,  two  important  advances  are  rapidly 
widening  the  scope  of  feasible  applications  of  the  LDA.  The  first  has  to  do  with  massive  paral¬ 
lelization  of  computer  hardware,  making  the  next  generation  of  computers  several  orders  of 
magnitude  more  powerful  than  present-day  ones.  Fortunately,  most  electronic  structure  methods 
are  quite  amenable  to  parallelization;  this  is  especially  so  with  respect  to  the  methods  described 
in  the  following  subsections.  The  second  has  to  do  with  significant  advances  in  algorithms  used 
to  solve  the  single-particle  Schrbdinger  equation. 
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6.1 


NEW  DEVELOPMENTS  IN  ELECTRONIC  STRUCTURE  ALGORITHMS 


Because  of  the  immense  potential  of  the  LDA  to  realistically  calculate  practical  matenal 
properties,  there  has  been  intensive  effort  to  develop  algorithms  to  solve  the  LDA  efficiently. 

Car  and  Parrinello^^  precipitated  the  present  efforts  by  proposing  an  algorithm  to  pursue  molecu¬ 
lar  dynamics  simulations  within  the  LDA,  by  introducing  a  technique  that  integrates  the  nuclear 
and  electronic  motion  into  a  single  step.  While  their  paper  attracted  a  great  deal  of  attention  for 
this  reason,  actually  the  most  important  result  of  their  paper  was  the  introduction  of  an  ingenious 
scheme  to  obtain,  for  a  plane-wave  (PW)  basis,  eigenvectors  of  the  Hamiltonian  much  more 
efficiently  than  previously.  They  exploited  the  fact  that  for  a  PW  basis,  the  kinetic-energy  pan  of 
the  Hamiltonian  is  diagonal  in  k-space  and  the  potential  is  diagonal  in  real-space.  Using  a  fast 
Fourier  transform  to  transform  a  trial  eigenvector  between  the  two  spaces,  they  could  operate  the 
Hamiltonian  on  a  trial  eigenvector  far  more  efficiently  than  before.  Using  an  iterative  scheme  to 
generate  the  eigenvectors  obviates  the  most  costly  of  the  n}  bottlenecks  in  solving  the  band- 
structure  problem.  To  iterate  M  trial  eigenvectors,  x,  in  a  PW  basis  of  dimension  N,  the  opera¬ 
tions  count  to  make  Hx  is  0(NM  In  AO.  In  this  procedure,  higher  lying  eigenstates  must  be 
orthogonalized  to  lower  ones  using  a  Gram-Schmidt  step,  which  requires  O(M^N)  operations.  M 
and  N  are  both  proportional  to  the  number  of  atoms  n,  so  the  bottleneck  remains;  on  a  PW 
basis,  however,  N»M  and  this  method  is  much  more  efficient  than  conventional  diagonalization 
schemes,  which  require  an  operations  count  of  order  0.  For  moderately  sized  systems,  the 
Gram-Schmidt  step  does  not  seem  to  dominate  the  total  time.  More  efficient  schemes  have  since 
been  introduced  that  allow  the  iteration  of  a  trial  eigenvector  to  proceed  to  the  eigenstate  in  fewer 

steps.28 


6.2  REAL-SPACE  HAMILTONIANS 

The  Car-Parrinello  method  (and  its  modern-day  descendants^^)  has  been  successful  for 
systems  of  moderate  size,  but  is  not  generally  applicable;  in  particular,  it  relies  on  pseudopoten¬ 
tial  approximations  that  make  it  ill  suited  to  many  applications  with  elements  outside  of  column,s 
HI- VII  in  the  Periodic  Table.  This  is  a  dynamic  and  rapidly  changing  subject,  to  which  we  can¬ 
not  do  justice  here.  We  will  briefly  describe  new  developments  we  have  been  pursuing,  that 
show  great  promise  for  yet  more  significant  improvements  in  efficiency.  The  barrier  is  not 
one  of  principle,  since  the  electronic  Hamiltonian  is  approximately  local:  well-separated  elec¬ 
trons  are  coupled  together  only  by  their  (classical)  electrostatic  interaction  and  it  is  unnecessary, 
in  principle,  to  couple  their  electronic  motion.  Thus  it  must  be  possible,  in  principle,  to  generate 
solutions  to  the  LDA  that  increase  only  linearly  with  the  number  of  atoms. 

By  construction,  the  PW  basis  is  nonlocal,  and  is  not  suited  to  techniques  that  exploit  the 
local  nature  of  the  Hamiltonian.  For  this,  a  real-space  basis  is  needed.  Real-space  basis  sets 
have  the  following  advantages:  (1)  They  are  tailored  to  the  potential  of  the  system;  thus  they  are 
much  smaller  and  more  efficient  than  a  PW  basis,  and  the  decomposition  of  eigenstates  into  local 
orbitals  helps  provide  physical  insight  into  the  problem.  (2)  Unlike  the  PW  basis,  they  are  gen¬ 
erally  applicable,  and  can  be  equally  well  adapted  to  “deep”  potentials,  such  as  are  encountered 
in  transition  metals.  (3)  A  real-space  approach  is  optimally  parallelizable.  As  is  well  known,  the 
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primar>'  obstacle  to  massively  paralic!  aiguruhnis  arises  iruin  delays  tn  inicrpt.Kti'.,*'  i  (‘wnum 
cations.  They  limit  the  overall  cfficictwy  ot  most  aigonihms  u*  a  small  Iracsion  uf  a  compuicr  s 
theoretical  efficiency.  Fortunately,  aigonihms  using  a  real- space  basis  can  tx-  tughiy  ^oaisc 
grained,  since  a  single  node  can  be  assigned  to  one  or  a  small  nurntscr  of  atoms,  and  requires 
minimal  communication  with  other  nodes  i4)  Using  a  real-space  basis,  it  is  fK'sssblc,  m  pnnci 
pie,  to  exploit  the  locality  of  the  Hamilioruan  and  make  the  total  computation  sirnc  increase  unis 
linearly  with  the  size  of  the  system.  Indeed,  the  recursion  method  is  a  strictly  real  spas c  tech 
nique,  for  which  the  compulation  time  increases  only  as  N  With  regard  to  parallch/juon,  the 
recursion  nr^thod  needs  no  interprocessor  communication  at  all  (with  the  minor  exception  of 
finding  the  Fermi  level  of  the  system  as  a  whole  at  the  end  of  an  itentiioni 

While  the  recursion  method  may  be  a  practicable  alternative  for  large  syssems,  finding  yes 
better  ways  to  diagonalize  a  real-space  Hamiltonian  is  one  key  pan  of  our  program  Ix>r  exam 
pie.  it  is  possible  to  adapt  some  of  the  conjugate-gradient  diagonalization  techniques*^  in  the  FV» 
basis  to  real-space  methods,  and  once  again  avoid  the  largest  of  the  bottlenecks  A  real  space 
Hamiltonian  is  local  and  thus  sparse  for  a  large  system  For  a  large  system,  the  operation  count 
to  make  Hx  intTcases  with  0{MN)  in  a  real  space  basis.  Both  a  FW  basis  and  a  real  space  basis 
have  a  Gram-Schmidt  orthogonalization  step  which  increases  as  0(.V/*.Vs  However  S  is  an  order 
of  magnitude  smaller  in  a  real-space  basis.  Still  other  approaches  can  be  considered,  and  study 
of  these  forms  one  key  aspect  of  our  research, 

6.3  CHOICE  OF  REAL-SPACE  HAMILTONIANS 

There  are  two  important  drawbacks  to  a  real-space  Hamiltoman.  The  first  is  that  actuaUy 
constructing  such  a  basis  |o}  and  calculating  matn.x  elements  for  it  is  in  general  quiic  difficult, 
indeed  much  more  complicated  than  for  a  PW  basis.  Any  mcthtxl  mu.st  earn'  out  the  following 
steps:  (1)  represent  the  charge  density  and  solve  Poisson’s  cuuation  for  the  clccirosiatic  poten 
tial;  (2)  evaluate  matrix  elements  of  the  Hamiltonian  j.  and  (.4)  sum  an  appropnaie  linear 

combination  of  wave  function  products  (0;  i(?j)  to  generate  the  output  density.  All  of  these  steps 
are  simply  done  using  plane  waves,  but  not  for  general  funebons  O.  A  second  major  draw  hack  iv 
that,  owing  to  the  minima!  basis,  intemuclear  forces  cannot  be  reliably  calculated  from  the 
Hellmann-Feynman  theorem,  even  though  knowledge  of  these  forces  is  essentia!  for  molecular 
dynamics  or  to  relax  a  large  configuration  of  atoms  to  their  equilibrium  positions 

One  way  to  tackle  these  obstacles  is  to  represent  the  charge  density  in  PW.  but  use  a  local 
basis  for  the  Hamiltonian.  Most  real-space  methods  use  this  technique,  and  we  employ  such  a 
procedure  in  the  present  work.  For  the  Hamiltonian,  we  employ  the  method  of  Linear  Muffin- 
Tin  Orbitals  invented  by  Andersen. Because  its  orbitals  are  tailored  to  the  potential  of  the 
solid,  the  LMTO  method  employs  the  most  intelligent  basis  set  of  all  electronic  structure  meth¬ 
ods.  Near  the  nuclei,  the  wave  function  is  represented  numerically  and  solves  the  Schrodinger 
equation  essentially  exactly.  In  the  interstitial  region,  the  kinetic  energy  of  the  envelope  funt 
tions  is  customarily  set  to  zero,  which  is  close  to  the  optimal  average  value  for  close-packed 
solids. 

An  alternative  technique,  very  recently  developed  by  Meihfessel  and  van  Schilfgaarde 
(manuscript  in  preparation)  uses  a  strictly  real-space  method,  in  which  both  the  basis  and  the 
density  are  represented  as  spatially  local,  atom-centered  functions.  This  new  technique  is  highly 
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efficient;  to  date  a  molecular  program  has  been  completed  in  which  it  is  possible  to  calculate  the 
eigenstates  of  small  molecules  in  a  few  seconds  on  a  workstation.  Here  we  sketch  only  the 
essentials  of  the  method;  the  paper  will  be  published  shonly.  In  any  real-space  method,  a  diffi¬ 
culty  immediately  ari.ses  when  attempting  to  represent  the  electron  density  p  which  is  calculated 
from  the  sum  of  products  of  basis  functions.  There  is  no  convenient  analytical  representation  of 
Hankel  function  products  that  the  LMTO  method  employs  (as  there  is  in  products  of  plane 
waves);  which  is  why  most  methods  resort  to  plane  waves  to  represent  the  density.  The  kernel  of 
the  present  method  centers  around  an  approximate — but  analytical — representation  of  the  prod¬ 
ucts  of  two  F  ike!  functions  centered  on  different  sites.  The  new  method  approximates  such  a 
product  as  a  linear  combination  of  Hankel  functions  centered  on  the  two  separate  sites.  (This 
approximation  is  done  once  and  for  all  by  tabulating  a  least-squares  fit  of  Hankel  function  prod¬ 
ucts,  so  that  mapping  a  product  of  two  functions  is  essentially  a  table  look-up.)  Thus  the  density 
p  can  be  mapped  into  a  linear  combination  of  atom-centered  Hankel  functions.  In  a  like  manner, 
matrix  elements  of  the  potential  are  easily  calculated,  once  again  using  the  mapping 

As  with  plane-wave  representations,  the  error  in  the  approximation  can  be  made  arbitrar¬ 
ily  small  by  increasing  the  number  of  functions  that  approximate  a  product.  The  key  advantage 
to  this  representation  is  that  atom-centered  Hankel  functions  are  strictly  local:  they  do  not 
extend  over  all  space  as  do  plane  waves.  Table  4  shows  some  calculated  results  for  selected 
small  molecules. 


Table  4 

COMPARISON  OF  STRUCTURAL  PROPERTIES  IN  SELECTED  MOLECULES 

(0  are  bond  angles  and  ware  vibrational  frequencies) 


<l(au) 

E(eV) 

e(deg) 

CO  (cnr^) 

Molecule 

Calc 

Expt 

Calc 

Expt 

Calc 

Expt 

Calc 

Expt 

N2 

2.074 

2.076 

11,34 

9.8 

— 

— 

2380 

2360 

Ha 

1.444 

1.410 

4.91 

4.5 

— 

— 

4270 

4400 

CU2 

4.124 

— 

2.97 

2.1 

— 

— 

291 

330 

H2O 

1.833 

1.811 

11.46 

9.5 

104.5 

104.5 

1560;  3700 

1600:3760 

H2S 

2.562 

2.543 

8.95 

7.5 

91.3 

92.3 

— 

— 

H2Se 

2.803 

2.779 

8.28 

— 

90.3 

91.0 

— 

— 

Also  with  this  method,  we  have  found  how  to  make  an  exact  differention  of  the  LDA 
Hamiltonian,  and  thus  obtain  a  proper  expression  for  the  forces  without  resorting  to  the 
Hellmann-Feynman  (force)  theorem.  The  usual  Hellmann-Feynman  theorem  uses  the  fact  that, 
at  self-consistency,  the  total  energy  is  stationary  to  first  order  with  respect  to  changes  in  the 
electron  density.  Minimal  basis  sets  incur  errors,  however,  because  the  basis  moves  with  virtual 
displacements  of  the  nuclei.  Our  force  theorem  overcomes  this  by  allowing  the  basis  and  charge 
density  to  shift  rigidly  with  nuclear  displacements.  We  can  do  so  because  they  are  represented  as 
a  superposition  of  atom-centered  functions.  Thus  we  obtain  additional  terms,  and  the  resulting 
expression  is  an  essentially  exact  expression  for  the  force. 
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To  summarize,  our  new  method  has  several  features  that  can  combine  the  best  of  many 
existing  electronic  structure  techniques.  It  uses  an  efficient,  minimal  basis,  but  at  the  same  time 
can  calculate  forces  precisely.  As  a  real-space  method,  it  can  be  exploited  to  render  it  highly 
efficient,  and  yet  accurate. 

As  yet,  we  only  have  a  real-space  molecule  program.  The  next  step  is  to  make  a  solid- 
state  program  based  on  the  same  principles.  We  anticipate  that  when  this  is  completed  in  the 
coming  months,  we  will  be  able  to  do  molecular  dynamics  simulations  of  complex  transition 
metal  and  metal-ceramic  systems  which  are  as  rapid  as  using  empirical  tight-binding,  but  are 
nonetheless  completely  first-principles  in  approach. 


I 
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